Study of the doubly charmed tetraquark \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{{\rm{T}}}}}}}_{{{{{{\rm{c}}}}}}{{{{{\rm{c}}}}}}}^{+}$$\end{document}Tcc+

Quantum chromodynamics, the theory of the strong force, describes interactions of coloured quarks and gluons and the formation of hadronic matter. Conventional hadronic matter consists of baryons and mesons made of three quarks and quark-antiquark pairs, respectively. Particles with an alternative quark content are known as exotic states. Here a study is reported of an exotic narrow state in the D0D0π+ mass spectrum just below the D*+D0 mass threshold produced in proton-proton collisions collected with the LHCb detector at the Large Hadron Collider. The state is consistent with the ground isoscalar \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{{\rm{T}}}}}}}_{{{{{{\rm{c}}}}}}{{{{{\rm{c}}}}}}}^{+}$$\end{document}Tcc+ tetraquark with a quark content of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{\rm{c}}}}}}{{{{{\rm{c}}}}}}\overline{{{{{{\rm{u}}}}}}}\overline{{{{{{\rm{d}}}}}}}$$\end{document}ccu¯d¯ and spin-parity quantum numbers JP = 1+. Study of the DD mass spectra disfavours interpretation of the resonance as the isovector state. The decay structure via intermediate off-shell D*+ mesons is consistent with the observed D0π+ mass distribution. To analyse the mass of the resonance and its coupling to the D*D system, a dedicated model is developed under the assumption of an isoscalar axial-vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{{\rm{T}}}}}}}_{{{{{{\rm{c}}}}}}{{{{{\rm{c}}}}}}}^{+}$$\end{document}Tcc+ state decaying to the D*D channel. Using this model, resonance parameters including the pole position, scattering length, effective range and compositeness are determined to reveal important information about the nature of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{{\rm{T}}}}}}}_{{{{{{\rm{c}}}}}}{{{{{\rm{c}}}}}}}^{+}$$\end{document}Tcc+ state. In addition, an unexpected dependence of the production rate on track multiplicity is observed.

H adrons with quark content other than that seen in mesons (q 1 q 2 ) and baryons (q 1 q 2 q 3 ) have been actively discussed since the birth of the quark model [1][2][3][4][5][6][7][8] . Since the discovery of the χ c1 (3872) state 9 many tetraquark and pentaquark candidates, listed in Table 1, have been observed [10][11][12][13][14][15][16][17][18][19] . For all but the X 0 (2900) and X 1 (2900) states the minimal quark content implies the presence of either a cc or bb quark-antiquark pair. The masses of many tetraand pentaquark states are close to mass thresholds, e.g., D ðÃÞ D ðÃÞ or B ðÃÞ B ðÃÞ , where D (*) or B (*) represents a hadron containing a charm or beauty quark, respectively. Therefore, these states are likely to be hadronic molecules 16,[20][21][22] where colour-singlet hadrons are bound by residual nuclear forces, such as the exchange of a pion or ρ meson 23 , similar to electromagnetic van der Waals forces attracting electrically neutral atoms and molecules. These states are expected to have a spatial extension significantly larger than a typical compact hadron. Conversely, the only hadron currently observed that contains a pair of c quarks is the Ξ þþ cc (ccu) baryon, a long-lived, weakly-decaying compact object 24,25 . The recently observed X(6900) structure in the J/ψJ/ψ mass spectrum 26 belongs to both categories simultaneously. Its proximity to the χ c0 χ c1 threshold could indicate a molecular structure 27,28 . Alternatively, it could be a compact object, where all four quarks are within one confinement volume and each quark interacts directly with the other three quarks via the strong force [29][30][31][32] .
The existence and properties of Q 1 Q 2 q 1 q 2 states with two heavy quarks and two light antiquarks have been widely discussed for a long time [33][34][35][36][37][38] . In the limit of large masses of the heavy quarks the corresponding ground state should be deeply bound. In this limit, the two heavy quarks Q 1 Q 2 form a point-like colourantitriplet object, analogous to an antiquark, and as a result, the Q 1 Q 2 q 1 q 2 system has similar degrees of freedom for its light quarks as an antibaryon with a single heavy quark, e.g., the Λ À c or Λ 0 b antibaryons. The beauty quark is considered heavy enough to sustain the existence of a bbud state that is stable with respect to the strong and electromagnetic interactions with a mass of about 200 MeV/c 2 below the BB * mass threshold. In the case of the bcud and ccud systems, there is currently no consensus in the literature whether such states exist and if their natural widths are narrow enough to allow for experimental observation. The theoretical predictions for the mass of the ccud ground state with spin-parity quantum numbers J P = 1 + and isospin I = 0, denoted hereafter as T þ cc , relative to the D *+ D 0 mass threshold lie in the range −300 < δm < 300 MeV/c 2 39-70 , where m D Ãþ and m D 0 denote the known masses of the D *+ and D 0 mesons 10 , with cd and cu quark content, respectively. The observation of a narrow state in the D 0 D 0 π + mass spectrum near the D *+ D 0 mass threshold, compatible with being a T þ cc tetraquark state with ccud quark content is reported in Ref. 71 .
In the work presented here, the properties of the T þ cc state are studied by constructing a dedicated amplitude model that accounts for the D *+ D 0 and D *0 D + decay channels. In addition, the mass spectra of other DD (*) and opposite-sign DD ðÃÞ combinations are explored. Furthermore, production-related observables, such as the event multiplicity and transverse momentum (p T ) spectra that are sensitive to the internal structure of the state, are discussed. This analysis is based on proton-proton (pp) collision data, corresponding to an integrated luminosity of 9 fb −1 , collected with the LHCb detector at centre-of-mass energies of 7, 8 and 13 TeV. The LHCb detector 72,73 is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks and is further described in Methods.

Results
T þ cc signal in the D 0 D 0 π + mass spectrum. The D 0 D 0 π + final state is reconstructed using the D 0 → K − π + decay channel with two D 0 mesons and a pion all produced promptly in the same pp collision. The inclusion of charge-conjugated processes is implied throughout the paper. The selection criteria are similar to those used in Refs. 74-77 and described in detail in Methods. The background not originating from true D 0 mesons is subtracted using an extended unbinned maximum-likelihood fit to the twodimensional distribution of the masses of the two D 0 candidates from selected D 0 D 0 π + combinations, see Methods and Supplementary Fig. 1a. The obtained D 0 D 0 π + mass distribution for selected D 0 D 0 π + combinations is shown in Fig. 1.
An extended unbinned maximum-likelihood fit to the D 0 D 0 π + mass distribution is performed using a model consisting of signal and background components. The signal component corresponds to the T þ cc ! D 0 D 0 π þ decay and is described as the convolution of the natural resonance profile with the detector mass resolution function. A relativistic P-wave two-body Breit-Wigner function F BW with a Blatt-Weisskopf form factor 78,79 is used in Ref. 71 as the natural resonance profile. That function, while sufficient to reveal the existence of the state, does not account for the resonance being in close vicinity of the D * D threshold. To assess the fundamental properties of resonances that are close to thresholds, advanced parametrisations ought to be used [80][81][82][83][84][85][86][87][88][89][90] . A unitarised Breit-Wigner profile F U , described in Methods Eq. (47), is used in this analysis. The function F U is built under two main assumptions.

Assumption 1
The newly observed state has quantum numbers J P = 1 + and isospin I = 0 in accordance with the theoretical expectation for the T þ cc ground state.

Assumption 2
The T þ cc state is strongly coupled to the D * D channel, which is motivated by the proximity of the T þ cc mass to the D * D mass threshold. Table 1 Tetra-and pentaquark candidates and their plausible valence quark content.
A study of the D 0 π + mass distribution for selected D 0 D 0 π + combinations in the region above the D *0 D + mass threshold and below 3.9 GeV/c 2 shows that approximately 90% of all D 0 D 0 π + combinations contain a true D *+ meson. Therefore, the background component is parameterised with a product of the twobody phase-space function Φ D Ãþ D 0 97 and a positive polynomial function P n , convolved with the detector resolution function R where n denotes the order of the polynomial function, n = 2 is used in the default fit. The D 0 D 0 π + mass spectrum with non-D 0 background subtracted is shown in Fig. 1 with the result of the fit using a model based on the F U signal profile overlaid. The fit gives a signal yield of 186 ± 24 and a mass parameter relative to the D *+ D 0 mass threshold, δm U of −359 ± 40 keV/c 2 . The statistical significances of the observed T þ cc ! D 0 D 0 π þ signal and for the δm U < 0 hypothesis are determined using Wilks' theorem to be 22 and 9 standard deviations, respectively.
The width of the resonance is determined by the coupling constant g for small values of g . With increasing g , the width increases to an asymptotic value determined by the width of the D *+ meson, see Methods and Supplementary Fig. 7 The mode relative to the D *+ D 0 mass threshold, δm, and the full width at half maximum (FWHM), w, for the F U profile are found to be δm ¼ À361 ± 40 keV=c 2 and w ¼ 47:8 ± 1:9 keV=c 2 , to be compared with those quantities determined for the F BW signal profile of δm ¼ À279 ± 59 keV=c 2 and w ¼ 409 ± 163 keV=c 2 . They appear to be rather different. Nonetheless, both functions properly describe the data given the limited sample size, and accounting for the detector resolution, and residual background. To quantify the impact of these experimental effects, two ensembles of pseudoexperiments are performed. Firstly, pseudodata samples are generated with a model based on the F U profile. The parameters used here are obtained from the default fit, and the size of the sample corresponds to the size of data sample. Each pseudodata sample is then analysed with a model based on the F BW function. The obtained mean and RMS values for the parameters δm BW and Γ BW over the ensemble are shown in Table 2. The mass parameter δm BW agrees well with the value determined from data 71 . The difference for the parameter Γ BW does not exceed one standard deviation. Secondly, an ensemble of pseudodata samples  Fig. 1 Distribution of D 0 D 0 π + mass. Distribution of D 0 D 0 π + mass where the contribution of the non-D 0 background has been statistically subtracted. The result of the fit described in the text is overlaid. Uncertainties on the data points are statistical only and represent one standard deviation, calculated as a sum in quadrature of the assigned weights from the background-subtraction procedure.   Table 2. These values agree well with the result of the default fit to data. The results of these pseudoexperiments explain the seeming inconsistency between the models and illustrate the importance of an accurate description of the detector resolution and residual background given the limited sample size.
Systematic uncertainties. Systematic uncertainties for the δm U parameter are summarised in Table 3 and described in greater detail below. The systematic uncertainty related to the fit model is studied using pseudoexperiments with a set of alternative parameterisations. For each alternative model, an ensemble of pseudoexperiments is performed with parameters obtained from a fit to data. A fit with the baseline model is performed on each pseudoexperiment, and the mean values of the parameters of interest are evaluated over the ensemble. The absolute values of the differences between these mean values and the corresponding parameter values obtained from the fit to data are used to assess the systematic uncertainty due to the choice of the fit model. The maximal value of such differences over the considered set of alternative models is taken as the corresponding systematic uncertainty. The following sources of systematic uncertainty related to the fit model are considered: • Imperfect knowledge of the detector resolution model. To estimate the associated systematic uncertainty a set of alternative resolution functions is tested: a symmetric variant of an Apollonios function 100 , a modified Gaussian function with symmetric power-law tails on both sides of the distribution 101,102 , a generalised symmetric Student's t-distribution 103,104 , a symmetric Johnson's S U distribution 105,106 , and a modified Novosibirsk function 107 .
• A small difference in the detector resolution between data and simulation. A correction factor of 1.05 is applied to account for known discrepancies in modelling the detector resolution in simulation. This factor was studied for different decays [94][95][96][108][109][110] and found to lie between 1.0 and 1.1. For decays with relatively low momentum tracks, this factor is close to 1.05, which is the nominal value used in this analysis. This factor is also cross-checked using large samples of D *+ → D 0 π + decays, where a value of 1.06 is obtained. To assess the systematic uncertainty related to this factor, detector resolution models with correction factors of 1.0 and 1.1 are studied as alternatives.
• Parameterisation of the background component. To assess the associated systematic uncertainty, the order of the positive polynomial function of Eq. (2) is varied. In addition, to estimate a possible effect from a small contribution from three-body D 0 D 0 π + combinations without an intermediate D *+ meson, a more general family of background models is tested where Φ D 0 D 0 π þ denotes the three-body phase-space function 111 . The functions B 0 , B 1 , B 3 and B 0 nm with n ≤ 2, m ≤ 1 are used as alternative models for the estimation of the systematic uncertainty.
• Values of the coupling constants for the D * → Dπ and D * → Dγ decays affecting the shape of the F U signal profile. These coupling constants are calculated from the known branching fractions of the D * → Dπ and D * → Dγ decays 10 , the measured natural width of the D *+ meson 10,112 and the derived value for the natural width of the D *0 meson 66,81,113 . To assess the associated systematic uncertainty, a set of alternative models built around the F U profiles, obtained with coupling constants varying within their calculated uncertainties, is studied.
• Unknown value of the g parameter. In the baseline fit the value of the g parameter is fixed to a large value. To assess the effect of this constraint the fit is repeated using the value of g ¼ 8:08 GeV, that corresponds to À2Δ log L ¼ 1 for the most conservative likelihood profile for g that accounts for the systematic uncertainty. The change of 7 keV/c 2 of the δm U parameter is assigned as the systematic uncertainty.
The calibration of the momentum scale of the tracking system is based upon large samples of B + → J/ψK + and J/ψ → μ + μ − decays 114 . The accuracy of the procedure has been checked using fully reconstructed B decays together with two-body ϒ(nS) and K 0 S decays and the largest deviation of the bias in the momentum scale of δα = 3 × 10 −4 is taken as the uncertainty 115 . This uncertainty is propagated for the parameters of interest using simulated samples, with momentum scale corrections of 1 ± δα ð Þ applied. Half of the difference between the obtained peak locations is taken as an estimate of the systematic uncertainty.
In the reconstruction the momenta of the charged tracks are corrected for energy loss in the detector material using the Bethe-Bloch formula 116,117 . The amount of the material traversed in the tracking system by a charged particle is known to have 10% accuracy 118 . To assess the corresponding uncertainty the magnitude of the calculated corrections is varied by ±10%. Half of the difference between the obtained peak locations is taken as an estimate of the systematic uncertainty due to energy loss corrections.
The mass of D 0 D 0 π + combinations is calculated with the mass of each D 0 meson constrained to the known value of the D 0 mass 10 . This procedure produces negligible uncertainties for the δm U parameter due to imprecise knowledge of the D 0 mass. However, the small uncertainty of 2 keV/c 2 for the known D *+ − D 0 mass difference 10,112,119 directly affects the values of these parameters and is assigned as a corresponding systematic uncertainty.
For the lower limit on the parameter g , only systematic uncertainties related to the fit model are considered. For each alternative model the likelihood profile curves are built and corresponding 90 and 95% CL lower limits are calculated using the procedure described above. The smallest of the resulting values is taken as the lower limit that accounts for the systematic uncertainty: g > 5:1 ð4:3Þ GeV at 90 (95%) CL.
Results. Studying the D 0 π + mass distribution for T þ cc ! D 0 D 0 π þ decays allows testing the hypothesis that the T þ cc ! D 0 D 0 π þ The total uncertainty is calculated as the sum in quadrature of all components. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30206-w decay proceeds through an intermediate off-shell D *+ meson. The background-subtracted D 0 π + mass distribution for selected D 0 D 0 π + candidates with the D 0 D 0 π + mass with respect to the D *+ D 0 mass threshold, δm D 0 D 0 π þ , below zero is shown in Fig. 3. Both D 0 π + combinations are included in this plot. The twodimensional distribution of the mass of one D 0 π + combination versus the mass of another D 0 π + combination is presented in Supplementary Fig. 10.
A fit is performed to this distribution with a model containing signal and background components. The signal component is derived from the A U amplitude, see Methods Eq. (48), and is convolved with a detector resolution for the D 0 π + mass. This detector resolution function is modelled with a modified Gaussian function with power-law tails on both sides of the distribution 101,102 and parameters taken from simulation. Similarly to the correction used for the D 0 D 0 π + mass resolution function R, the width of the Gaussian function is corrected by a factor of 1.06 which is determined by studying large samples of D *+ → D 0 π + decays. The RMS of the resolution function is around 220 keV/c 2 . The shape of the background component is derived from data for δm D 0 D 0 π þ >0:6 MeV=c 2 . The fit results are overlaid in Fig. 3. The background component vanishes in the fit, and the D 0 π + spectrum is consistent with the hypothesis that the T þ cc ! D 0 D 0 π þ decay proceeds through an intermediate off-shell D *+ meson. This in turn favours the 1 + assignment for the spinparity of the state. Due to the proximity of the observed T þ cc signal to the D *+ D 0 mass threshold, and the small energy release in the D *+ → D 0 π + decay, the D 0 D 0 mass distribution from theT þ cc ! D 0 D 0 π þ decay forms a narrow peak just above the D 0 D 0 mass threshold. In a similar way, a peaking structure in the D + D 0 mass spectrum just above the D + D 0 mass threshold is expected from T þ cc ! D þ D 0 π 0 and T þ cc ! D þ D 0 γ decays, both proceeding via off-shell intermediate D *+ D 0 and D *0 D + states. The D 0 D 0 and D + D 0 final states are reconstructed and selected similarly to the D 0 D 0 π + final state, where the D + → K − π + π − decay channel is used. The background-subtracted D 0 D 0 and D + D 0 mass distributions are shown in Fig. 4 (top), where narrow structures are clearly visible just above the DD thresholds. Fits to these distributions are performed using models consisting of two components: a signal component F DD described in Methods Eqs. (70) and (71) and obtained via integration of the matrix elements for the T þ cc ! DDπ=γ decays with the F U profile, and a background component, parameterised as a product of the two-body phasespace function Φ DD and a positive linear function P 1 . The fit results are overlaid in Fig. 4 (top). The signal yields in the D 0 D 0 and D + D 0 spectra are found to be 263 ± 23 and 171 ± 26, respectively. The statistical significance of the observed T þ cc ! D 0 D 0 X and T þ cc ! D þ D 0 X signals, where X stands for nonreconstructed pions or photons, is estimated using Wilks' theorem 120 and is found to be in excess of 20 and 10 standard deviations, respectively. The relative yields for the signals observed in the D 0 D 0 π + , D 0 D 0 and D 0 D + mass spectra agree with the expectations of the model described in Methods where the decay of an isoscalar T þ cc state via the D * D channel with an intermediate off-shell D * meson is assumed.
The observation of the near-threshold signals in the D 0 D 0 and D + D 0 mass spectra, along with the signal shapes and yields, all agree with the isoscalar T þ cc hypothesis for the narrow signal observed in the D 0 D 0 π + mass spectrum. However, an alternative interpretation could be that this state is the I 3 = 0 component of â T cc isotriplet (T 0 cc ,T þ cc ,T þþ cc ) with ccuu, ccud and ccdd quark content, respectively. Assuming that the observed peak corresponds to theT þ cc component and using the estimates for theT cc mass splitting from Methods Eqs. (85) and (86), the masses of thê T 0 cc andT þþ cc states are estimated to be slightly below the D 0 D *0 and slightly above the D + D *+ mass thresholds, respectively: With these mass assignments, assuming equal production of all threeT cc components, theT 0 cc state would be an extra narrow state that decays into the D 0 D 0 π 0 and D 0 D 0 γ final states via an off-shell D *0 meson. These decays would contribute to the narrow near-threshold enhancement in the D 0 D 0 spectrum, and increase the signal in the D 0 D 0 mass spectrum by almost a factor of three. TheT þþ cc state would decay via an on-shell D *+ meson T þþ cc ! D þ D Ãþ ; therefore, it could be a relatively wide state, with a width up to a few Me 121 . Therefore, it would manifest itself as a peak with a moderate width in the D + D 0 π + mass spectrum with a yield comparable to that of theT þ cc ! D 0 D 0 π þ decays. In addition, it would contribute to the D + D 0 mass spectrum, tripling the contribution from theT þ cc decays. However, due to the larger mass of theT þþ cc state and its larger width, this contribution should be wider, making it more difficult to disentangle from the background. Finally, theT þþ cc state would make a contribution to the D + D + spectrum with a yield similar to the contribution from T þ cc ! D 0 D þ π 0 =γ decays to the D 0 D + spectrum, but wider. The mass spectra for D + D + and D + D 0 π + combinations are shown in Fig. 4 (bottom).
Neither distribution exhibits any narrow signal-like structure. Fits to these spectra are performed using the following background-only functions: The results of these fits are overlaid in Fig. 4 (bottom). The absence of any signals in the D + D + and D + D 0 π + mass spectra is therefore a strong argument in favour of the isoscalar nature of the observed peak in the D 0 D 0 π + mass spectrum. The interference between two virtual channels for the T þ cc ! D 0 D 0 π þ decay, corresponding to two amplitude terms, see Methods Eq. (35), is studied by setting the term proportional to C in Methods Eq. (39) to be equal to zero. This causes a 43% reduction in the decay rate, pointing to a large interference. The same procedure applied to the T þ cc ! D þ D 0 π 0 decays gives the contribution of 45% for the interference between the D Ãþ ! D þ π 0 À Á D 0 and D Ã0 ! D 0 π 0 À Á D þ channels. For T þ cc ! D þ D 0 γ decays the role of the interference between the D Ãþ ! D þ γ À Á D 0 and D Ã0 ! D 0 γ À Á D þ channels is estimated by equating to zero the F þ F Ã 0 and F Ã þ F 0 terms in Methods Eqs. (45) and (46). The interference contribution is found to be 33%.
Using the model described earlier and results of the fit to the D 0 D 0 π + mass spectrum, the position of the amplitude poleŝ in the complex plane, responsible for the appearance of the narrow structure in the D 0 D 0 π + mass spectrum is determined. The pole parameters, mass m pole and width Γ pole , are defined through the pole locationŝ as ffif The pole locationŝ is a solution to the equation where A I I U ðsÞ denotes the amplitude on the second Riemann sheet defined in Methods Eq. (64). For large coupling g the position of the resonance pole is uniquely determined by the parameter δm U , i.e., the binding energy and the width of the D *+ meson. Figure 5 shows the complex plane of the δ ffiffi s p variable, defined as All possible positions of the pole for g ) m D 0 þ m D Ãþ are located on a red dashed curve in Fig. 5. The behaviour of the curve can be understood as follows: with an increase of the binding energy (distance to the D *+ D 0 mass threshold), the width gets narrower; and when the parameter δm U approaches zero, the pole touches the D 0 D *+ cut and moves to the other complex sheet, i.e., the state becomes virtual. For smaller values of g , the corresponding to the openings of the D 0 D + γ, D 0 D + π 0 and D 0 D 0 π + decay channels, are outside of the displayed region.
pole is located between the limiting curve and the ℜ s = 0 line. The pole parameters are found to be where the first uncertainty is due to the δm U parameter and the second is due to the unknown value of the g parameter. The peak is well separated from the D *+ D 0 threshold in the D 0 D 0 π + mass spectrum. Hence, as for an isolated narrow resonance, the parameters of the pole are similar to the visible peak parameters, namely the mode δm and FWHM w.
The systematic uncertainties quoted here do not account for the possibility that any of the underlying assumptions on which the model is built are not valid. For example, as shown earlier the data are consistent with a wide range of FHWM w values for the signal profile. Therefore the pole width Γ pole is based mainly on the T þ cc amplitude model and the value of the m U parameter determined from the fit to the D 0 D 0 π + mass spectrum.
A study of the behaviour of the A U ðsÞ amplitude in the vicinity of the D *+ D 0 mass threshold leads to the determination of the low-energy scattering parameters, namely the scattering length, a, and the effective range, r. These parameters are defined via the coefficients of the first two terms of the Taylor expansion of the inverse non-relativistic amplitude 122 , i.e., where k is the wave number.
Typically, a non-vanishing imaginary part of the scattering length indicates the presence of inelastic channels 123 ; however, in this case the non-zero imaginary part is related to the lower threshold, T þ cc ! D 0 D 0 π þ , and is determined by the width of the D *+ meson. The real part of the scattering length a is negative indicating attraction. This can be interpreted as the characteristic size of the state 16 , R a À< a ¼ 7:16 ± 0:51 fm: For the A U amplitude the effective range r is non-positive and proportional to g À2 , see Methods Eq. (76). Its value is consistent with zero for the baseline fit. An upper limit on the −r value is set as 0 ≤ À r < 11:9 ð16:9Þ fm at 90 ð95Þ% CL: The Weinberg compositeness criterion 124,125 makes use of the relation between the scattering length and the effective range to construct the compositeness variable Z, Another estimate of the characteristic size is obtained from the value of the binding energy ΔE. Within the interpretation of the T þ cc state as a bound D *+ D 0 molecular-like state, the binding energy is ΔE = − δm U . The characteristic momentum scale γ 16 is estimated to be where μ is the reduced mass of the D *+ D 0 system. This value of the momentum scale in turn corresponds to a characteristic size R ΔE of the molecular-like state, which is consistent with the R a estimate from the scattering length.
For high-energy hadroproduction of a state with such a large size, R a or R ΔE , one expects a strong dependency of the production rate on event multiplicity, similar to that observed for the χ c1 (3872) state 126  production cross-section is suppressed with respect to the conventional charmonium state ψ(2S) at large track multiplicities 126 . It is noteworthy that the track multiplicity distribution for the T þ cc state differs from that of the low-mass D 0 D 0 pairs, in particular, no suppression at large multiplicity is observed. A p value for the consistency of the track multiplicity distributions for T þ cc production and low-mass D 0 D 0 pairs is found to be 0.1%. It is interesting to note that the multiplicity distribution for T þ cc production and the one for D 0 D 0 -pairs with 3:75 < m D 0 D 0 < 3:87 GeV=c 2 are consistent with a corresponding p value of 12%. The similarity between T þ cc production, which is inherently a single parton scattering process, and the distribution for process dominated by a double-parton scattering is surprising.
The transverse momentum spectrum for the T þ cc state is compared with those for the low-mass D 0 D 0 and D 0 D 0 pairs in Fig. 6 (right). The p values for the consistency of the p T spectra for the T þ cc state and low-mass D 0 D 0 pairs are 1.4%, and 0.02% for low-mass D 0 D 0 pairs. More data are needed for further conclusions.
The background-subtracted D 0 D 0 mass distribution in a wider mass range is shown in Fig. 7 together with a similar distribution for D 0 D 0 pairs. In the D 0 D 0 mass spectrum the near-threshold enhancement is due to χ c1 ð3872Þ ! D 0 D 0 π 0 and χ c1 ð3872Þ ! D 0 D 0 γ decays via intermediate D *0 mesons 77 . This structure is significantly wider than the structure in the D 0 D 0 mass spectrum from T þ cc ! D 0 D 0 π þ decays primarily due to the larger natural width and smaller binding energy for the χ c1 (3872) state 94,95 . With more data, and with a better understanding of the dynamics of χ c1 ð3872Þ ! D 0 D 0 π 0 =γ decays, and therefore of the corresponding shape in the D 0 D 0 mass spectrum, it will be possible to estimate the relative production rates for the T þ cc and χ c1 (3872) states. Background-subtracted D 0 D 0 π + and D 0 D + mass distributions together with those for D 0 D 0 π þ and D 0 D − are shown in Supplementary Figs. 3 and 4.

Discussion
The exotic narrow tetraquark state T þ cc observed in Ref. 71 is studied using a dataset corresponding to an integrated luminosity of 9 fb −1 , collected by the LHCb experiment in pp collisions at centre-ofmass energies of 7, 8 and 13 TeV. The observed D 0 π + mass distribution indicates that the T þ cc ! D 0 D 0 π þ decay proceeds via an intermediate off-shell D *+ meson. Together with the proximity of the state to the D * D 0 mass threshold, this favours the spin-parity quantum numbers J P to be 1 + . Narrow near-threshold structures are observed in the D 0 D 0 and D 0 D + mass spectra with high significance. These are found to be consistent with originating from off-shell T þ cc ! D Ã D decays followed by the D * → Dπ and D * → Dγ decays. No signal is observed in the D + D 0 π + mass spectrum, and no structure is observed in the D + D + mass spectrum. These non-observations provide a strong argument in favour of the isoscalar nature for the observed state, supporting its interpretation as the isoscalar J P = 1 + ccud-tetraquark ground state. A dedicated unitarised three-body Breit-Wigner amplitude is built on the assumption of strong isoscalar coupling of the axial-vector T þ cc state to the D * D channel. This assumption is supported by the data; however, alternative models are not excluded by the distributions studied in this analysis. Probing alternative models and the validity of the underlying assumptions of this analysis will be a subject for future studies.
Using the developed amplitude model, the mass of the T þ cc state, relative to the D *+ D 0 mass threshold, is determined to be where the first uncertainty is statistic and the second systematic. The lower limit on the absolute value of the coupling constant of the T þ cc state to the D * D system is g > 5:1 ð4:3Þ GeV at 90 ð95Þ % CL: Using the same model, the estimates for the scattering length a, effective range r, and the compositeness, Z are obtained from the low-energy limit of the amplitude to be The characteristic size calculated from the binding energy is R ΔE = 7.49 ± 0.42 fm. This value is consistent with the estimation from the scattering length, R a = 7.16 ± 0.51 fm. Both R ΔE and R a correspond to a spatial extension significantly exceeding the typical scale for heavy-flavour hadrons. Within this model the resonance pole is found to be located on the second Riemann sheet with respect to the D 0 D 0 π + threshold, atŝ ¼ m pole À i 2 Γ pole , where  δm pole ¼ À360 ± 40 þ4 À0 keV=c 2 ; ð26Þ where the first uncertainty accounts for statistical and systematic uncertainties for the δm U parameters, and the second is due to the unknown value of the g parameter. The pole position, scattering length, effective range and compositeness form a complete set of observables related to the T þ cc ! D 0 D 0 π þ reaction amplitude, which are crucial for inferring the nature of the T þ cc tetraquark. Unlike in the prompt production of the χ c1 (3872) state, no suppression of the T þ cc production at high track multiplicities is observed relative to the low-mass D 0 D 0 pairs. The observed similarity with the multiplicity distribution for the low-mass D 0 D 0 production process, that is presumably double-partonscattering dominated, is unexpected. In the future with a larger dataset and including other decay modes, e.g., D 0 → K − π + π + π − , detailed studies of the properties of this new state and its production mechanisms could be possible.
In conclusion, the T þ cc tetraquark observed in D 0 D 0 π + decays is studied in detail, using a unitarised model that accounts for the relevant thresholds by taking into account the D 0 D 0 π + and D 0 D + π 0 (γ) decay channels with intermediate D * resonances. This model is found to give an excellent description of the D 0 π + mass distribution in the T þ cc ! D 0 D 0 π þ decay and of the threshold enhancements observed in the D 0 D 0 and D 0 D + spectra. Together with the absence of a signal in the D 0 D + and D + D 0 π + mass distributions this provides a strong argument for interpreting the observed state as the isoscalar T þ cc tetraquark with spin-parity J P = 1 + . The precise T þ cc mass measurement will rule out or improve on a considerable range of theoretical models on heavy quark systems. The determined pole position and physical quantities derived from low-energy scattering parameters reveal important information about the nature of the T þ cc tetraquark. In addition, the counterintuitive dependence of the production rate on track multiplicity will pose a challenge for theoretical explanations.

Methods
Experimental setup. The LHCb detector 72,73 is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary pp collision vertex, the impact parameter, is measured with a resolution of (15 + 29/p T ) μm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors 128 . Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction.
Simulation. Simulation is required to model the effects of the detector acceptance, resolution, and the efficiency of the imposed selection requirements. In the simulation, pp collisions are generated using PYTHIA 129 with a specific LHCb configuration 130 . Decays of unstable particles are described by EVTGEN 131 , in which final-state radiation is generated using PHOTOS 132 . The interaction of the generated particles with the detector, and its response, are implemented using the GEANT4 toolkit 133,134 as described in Ref. 135 .
Event selection. The D 0 D 0 , D 0 D + and D 0 D 0 π + final states are reconstructed using the D 0 → K − π + and D + → K − π + π + decay channels. The selection criteria are similar to those used in Refs. [74][75][76][77] . Kaons and pions are selected from well-reconstructed tracks within the acceptance of the spectrometer that are identified using information from the ring-imaging Cherenkov detectors. The kaon and pion candidates that have transverse momenta larger than 250 MeV/c and are inconsistent with being produced at a pp interaction vertex are combined together to form D 0 and D + candidates, referred to as D hereafter. The resulting D candidates are required to have good vertex quality, mass within ±65 and ±50 MeV/c 2 of the known D 0 and D + masses 10 , respectively, transverse momentum larger than 1 GeV/c, decay time larger than 100 μm/c and a momentum direction that is consistent with the vector from the primary to secondary vertex. Selected D 0 and D + candidates consistent with originating from a common primary vertex are combined to form D 0 D 0 and D 0 D + candidates. The resulting D 0 D 0 candidates are combined with a pion to form D 0 D 0 π + candidates. At least one of the two D 0 π + combinations is required to have good vertex quality and mass not exceeding the known D *+ mass by more than 155 MeV/c 2 . For each D 0 D 0 , D 0 D + and D 0 D 0 π + candidate a kinematic fit 136 is performed. This fit constrains the mass of the D candidates to their known values and requires both D mesons, and a pion in the case of D 0 D 0 π + , to originate from the same primary vertex. A requirement is applied to the quality of this fit to further suppress combinatorial background and reduce background from D candidates produced in two independent pp interactions or in the decays of beauty hadrons 74 . To suppress background from kaon and pion candidates reconstructed from a common track, all track pairs of the same charge are required to have an opening angle inconsistent with zero and the mass of the combination must be inconsistent with the sum of the masses of the two constituents. For cross-checks additional final states D + D + , D + D 0 π + , D 0 D 0 , D 0 D − and D 0 D 0 π þ are reconstructed, selected and treated in the same way.
Non-D background subtraction. Two-dimensional distributions of the mass of one D candidate versus the mass of the other D candidate from selected D 0 D 0 π + , D 0 D 0 and D 0 D + combinations are shown in Supplementary Fig. 1. These distributions illustrate the relatively small combinatorial background levels due to fake D candidates. This background is subtracted using the sPlot technique 137 , which is based on an extended unbinned maximum-likelihood fit to these two-dimensional distributions with the function described in Ref. 74 . This function consists of four components: • a component corresponding to genuine D 1 D 2 pairs and described as a product of two signal functions, each parameterised with a modified Novosibirsk function 107 ; • two components corresponding to combinations of one of the D mesons with combinatorial background, described as a product of the signal function and a background function, which is parameterised with a product of an exponential function and a positive first-order polynomial; • a component corresponding to pure background pairs and described by a product of exponential functions and a positive two-dimensional nonfactorisable second-order polynomial function.
Based on the results of the fit, each candidate is assigned a positive weight for being signal-like or a negative weight for being background-like, with the masses of the two D 0 candidates as discriminating variables. The D 0 D 0 π + mass distributions for each of the subtracted background components are presented in Supplementary  Fig. 2. where fit results with background-only functions B 0 10 , defined in Eq. (3) are overlaid.
Resolution model for the D 0 D 0 π + mass. In the vicinity of the D *+ D 0 mass threshold the resolution function R for the D 0 D 0 π + mass is parametrised with the sum of two Gaussian functions with a common mean. The widths of the Gaussian functions are σ 1 = 1.05 × 263 keV/c 2 and σ 2 = 2.413 × σ 1 for the narrow and wide components, respectively, and the fraction of the narrow Gaussian is α = 0.778. The parameters α and σ 1,2 are taken from simulation, and σ 1,2 are corrected with a factor of 1.05 that accounts for a small difference between simulation and data for the mass resolution [94][95][96] . The RMS of the resolution function is around 400 keV/c 2 .
Matrix elements for T þ cc ! DDπ=γ decays. Assuming isospin symmetry, the isoscalar vector state T þ cc that decays into the D * D final state can be expressed as Therefore, the S-wave amplitudes for the T þ cc ! D Ãþ D 0 and T þ cc ! D Ã0 D þ decays have different signs where g is a coupling constant, ϵ T þ cc is the polarisation vector of the T þ cc particle and ϵ D Ã is the polarisation vector of D * meson, and the upper and lower Greek indices imply the summation in the Einstein notation. The S-wave (corresponding to orbital angular momentum equal to zero) approximation is valid for a nearthreshold peak. For T þ cc masses significantly above the D * D threshold, higher-order waves also need to be considered. The amplitudes for the D * → Dπ decays are written as where f denotes a coupling constant, and p D stands for the momentum of the D meson. The amplitude for the D * → Dγ decays is where h denotes a coupling constant, μ stands for the magnetic moment for D * → Dγ transitions, p D Ã and p γ are the D * -meson and photon momenta, respectively, and ϵ γ is the polarisation vector of the photon. The three amplitudes for T þ cc ! πDD and T þ cc ! γDD decays are where s ij ¼ p 2 ij ¼ ðp i þ p j Þ 2 and the F functions that denote the Breit-Wigner amplitude for the D * mesons are A small possible distortion of the Breit-Wigner shape of the D * meson due to three-body final-state interactions is neglected in the model. The impact of the energy-dependence of the D * meson self-energy is found to be insignificant. The decays of the T þ cc state into the D + D + π − final state via off-shell D *0 → D + π − decays are highly suppressed and are not considered here. The last terms in Eqs. (35) and (36) imply the same amplitudes with swapped momenta.
The T þ cc state is assumed to be produced unpolarised; therefore, the squared absolute value of the decay amplitudes with pions in the final state, averaged over the initial spin state are and λ x; y; z À Á stands for the Källén function 97 . The additional factor of 2! in the denominator of Eq. (39) is due to the presence of two identical particles (D 0 ) in the final state. The squared absolute values of the decay amplitude with a photon in the final state, averaged over the initial spin state is The coupling constants f and h for the D * → Dπ and D * → Dγ decays are calculated using Eqs. The magnetic moment μ + is taken to be 1 and the ratio of magnetic moments μ 0 / μ + is calculated according to Refs. [138][139][140] .
Unitarised Breit-Wigner shape. A unitarised three-body Breit-Wigner function is defined as where f 2 D 0 D 0 π þ ; D 0 D þ π 0 ; D 0 D þ γ È É denotes the final state. The decay matrix element for each channel integrated over the three-body phase-space is denoted by where M For large values of s, in excess of s * , such as ffiffiffiffi where Φ D Ã D ðsÞ denotes the two-body phase-space function, the constants c 1 ,c 2 and c 3 are chosen to ensure the continuity of the functions ϱ f (s), and a value of ffiffiffiffi s Ã p ¼ 3:9 GeV=c 2 is used. The functions ϱ f (s) are shown in Supplementary Fig. 5. The complex-valued widthΓðsÞ is defined via the self-energy function Σ(s) 141 where g 2 is again factored out for convenience. The imaginary part of Σ(s) for real physical values of s is computed through the optical theorem as half of the sum of the decay probability to all available channels 142 : The real part of the self-energy function is computed using Kramers-Kronig dispersion relations with a single subtraction 143,144 , ξðsÞ ¼ s 2π p:v:  Supplementary Fig. 6.
Alternatively, the isoscalar amplitude A U is constructed using the K-matrix approach 145 with two coupled channels, D *+ D 0 and D *0 D + . The relation reads: where a production vector P and an isoscalar potential K are defined as The propagation matrix G describes the D * D → D * D rescattering via the virtual loops including the one-particle exchange process 91 and expressed in a symbolical way in Supplementary Eq. (1). where suppressed D *0 → D + π − transition is neglected. The D * + D 0 ↔ D *0 D + rescattering occurs due to non-diagonal element of the K-matrix (contact interaction) and non-diagonal elements of the G matrix (long-range interaction). The matrix G and the self-energy function Σ(s) from Eqs. (54) and (56), are related as Similar to the Flatté function 98 for large values of the g parameter, the F U signal profile exhibits a scaling property 94,99 . For large values of the g parameter the width approaches asymptotic behaviour, see Supplementary Fig. 7. The unitarised threebody Breit-Wigner function F U for T þ cc ! D 0 D 0 π þ decays with parameters m U and g obtained from the fit to data is shown in Supplementary Fig. 8. The inset illustrates the similarity of the profile with the single-pole profile in the vicinity of the pole ffif where m and w are the mode and FWHM, respectively.
Analytic continuation. Equation (56) defines Σ(s) and the amplitude A U ðsÞ for real values of s. Analytic continuation to the whole first Riemann sheet is calculated as where the integral is understood as in Eq. (58). The search for the resonance pole requires knowledge of the amplitude on the second Riemann sheet denoted by A I I U . According to the optical theorem 142 , the discontinuity of the inverse amplitude across the unitarity cut is given by i g 2 ϱ tot ðsÞ: For the complex-values s, the analytic continuation of ϱ tot (s) is needed: the phase-space integral in Eq. (49) is performed over a two-dimensional complex manifold D (see discussion on the continuation in Ref. 146 ): The integration is performed along straight lines connecting the end points in the complex plane.
2. A power-law cutoff f P C ðsÞ defined as Fits to the background-subtracted D 0 D 0 π + mass spectrum using a signal profile of the form F U ðsÞ f C ðsÞ show that the parameter δm U is insensitive to the choice of cutoff function when x c ≥ m D Ã0 þ m D þ and σ c ≥1 MeV/c 2 . The power-law cutoff function f P C ðsÞ with parameters x c ¼ m D Ã0 þ m D þ and σ c = 1 MeV/c 2 is chosen. The shapes for the D 0 D 0 and D + D 0 mass distributions are defined as Low-energy scattering amplitude. The unitarized Breit-Wigner amplitude is formally similar to the low-energy expansion given by Eq. (13) once the factor 1 2 g 2 is divided out 2 The function iϱ tot (s) matches ik up to a slowly varying energy factor that can be approximated by a constant in the threshold region. The proportionality factor w has the dimension of an inverse mass and is found by matching the decay probability to the two-body phase-space expression: where c 1 is a coefficient computed in Eq. (50). The comparison of A À1 NR and A À1 U 2w= g 2 that validates the matching is shown in Supplementary Fig. 9.
The inverse scattering length is defined as the value of the amplitude in Eq. (72) at the D *+ D 0 threshold: The imaginary part is fully determined by the available decay channels, while the real part depends on the constant ξðm 2 U Þ adjusted in the fit. The quadratic term, k 2 in Eq. (72), corresponds to the linear correction ins since k 2 = (s − s th )/4 for the non-relativistic case. Hence, the slope of the linear term in the A À1 U amplitude is related to the effective range as follows: Mass splitting for theT cc isotriplet. While the degrees of freedom of the light diquark for the isoscalar T þ cc state are similar to those for the Λ À c state, for theT cc isotriplet (T 0 cc ,T þ cc ,T þþ cc ) the light diquark degrees of freedom would be similar to those for the Σ c (anti)triplet. Assuming that the difference in the light quark masses, the Coulomb interaction of light quarks in the diquark, and the Coulomb interaction of the light diquark with the c-quark are responsible for the observed mass splitting in the Σ c isotriplet, the masses for the Σ c states can be written as where m Σ is a common mass parameter; the second and third terms describe the contribution from the light quark masses, m u and m d , into the mass splitting; terms proportional to a describe Coulomb interactions of light quarks in the diquark; terms proportional to b describe the Coulomb interactions of the diquark with the c-quark; and q q denotes the charge of the q-quark. Similar expressions can be written for theT cc isotriplet: where mT cc is the common mass parameter, q q ¼ Àq q and q cc = 2q c is the charge of a cc diquark. Using the known masses of the light quarks and Σ c states 10 and taking a 0 ¼ a and b 0 ¼ b, the mass splitting for theT cc isotriplet is estimated to be The validity of this approach is tested by comparing the calculated mass splitting between Σbp and Σbm states of −6.7 ± 0.7 MeV/c 2 with the measured value of −5.1 ± 0.2 MeV/c 210 . Based on the small observed difference, an additional uncertainty of 0.8 MeV/c 2 is added in quadrature to the results from Eqs. (83) and (84), and finally one gets These results agree with the assigned uncertainty with results based on a more advanced model from Ref. 147 .

Data availability
LHCb data used in this analysis will be released according to the LHCb external data access policy, which can be downloaded from https://opendata.cern.ch/record/410/files/ LHCb-Data-Policy.pdf. The raw data in all of the figures of this manuscript can be downloaded from https://cds.cern.ch/record/2780001, where no access codes are required. In addition, the unbinned background-subtracted data, shown in Figs. 1, 3 and 4 have been added to the HEPDATA record at https://www.hepdata.net/record/ins1915358.

Code availability
LHCb software used to process the data analysed in this manuscript is available at GITLAB repository. The specific software used in data analysis is available at ZENODO  166. CDF Collaboration, Aaltonen, T. et al. Evidence for a narrow near-threshold structure in the J/ψ ϕ mass spectrum in B + → J/ψ ϕK + decays.